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Abstract 

In  this  paper  we  explicitly  construct  the  entropy  solutions  for  the  Lighthill-Whitham- 
Richards  (LWR)  traffic  flow  model  with  a  flow-density  relationship  q{p)  which  is  piece- 
wise  quadratic,  continuous,  concave,  but  not  differentiable  at  the  junction  points  where 
two  quadratic  polynomials  meet,  and  with  piecewise  linear  initial  condition  and  piecewise 
constant  boundary  conditions.  As  observed  traffic  flow  data  can  be  well  fitted  with  such 
continuous  piecewise  quadratic  functions,  the  explicitly  constructed  solutions  provide  a  fast 
and  accurate  solution  tool  which  may  be  used  for  predicting  traffic  or  as  a  diagnosing  tool  to 
test  the  performance  of  numerical  schemes.  We  implement  these  explicit  entropy  solutions 
for  three  representative  traffic  flow  cases  and  also  compare  them  with  numerical  solutions 
obtained  by  a  high  order  weighted  essentially  non-oscillatory  (WENO)  scheme. 
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1  Introduction 


One  of  the  earliest  studies  of  traffic  flow  theory  began  when  Greenshields  (1934)  measured 
the  traffic  flow  and  speed  on  a  highway  and  identified  a  linear  relationship  between  speed  and 
density.  To  describe  the  dynamic  characteristics  of  traffic  on  a  homogeneous  and  unidirec¬ 
tional  highway,  Lighthill  and  Whitham  (1955)  and  Richards  (1956)  independently  proposed 
a  macroscopic  model  of  traffic  flow,  which  is  now  known  as  the  LWR  model  in  the  litera¬ 
ture  of  traffic  flow  theory.  Since  then,  a  substantial  amount  of  work  has  been  conducted  to 
improve  the  modeling  approach  of  traffic  flows,  which  can  be  broadly  classified  as  hydro- 
dynamic  models  (Daganzo,  1994;  Kerner  and  Konhauser,  1994;  Newell,  1993;  Payne,  1979; 
Zhang,  1998;  Wong  and  Wong,  2002a),  gas  kinetic  models  (Helbing,  1996;  Hoogendoorn 
and  Bovy,  2000;  Nelson  and  Sopasakis,  1998;  Prigogine  and  Herman,  1971;  Pavieri-Fontana, 
1975),  and  cellular  automation  models  (Biham  et  ah,  1992;  Cuesta  et  al.,  1993;  Krauss  et 
ah,  1997;  Nagatani,  1993;  Nagel  and  Schreckenberg,  1992).  Although  advances  have  been 
made  in  many  directions,  the  LWR  model  is  still  widely  used  for  the  modeling  of  traffic  flow, 
because  of  its  simplicity  and  good  explanatory  power  to  understand  the  qualitative  behavior 
of  road  traffic.  The  results  that  are  obtained  from  the  LWR  model  are  generally  adequate 
for  many  applications  such  as  traffic  management  and  control  problems. 

The  LWR  model  is  formulated  as  a  scalar  hyperbolic  conservation  law  and  is  often  solved 
by  finite  difference  methods  (Daganzo,  1995;  LeVeque,  1992;  Lebacque,  1996;  Michalopoulos 
et  ah,  1984).  The  main  difficulty  in  designing  efficient  and  high  order  finite  difference  methods 
for  the  LWR  model  or  in  general  for  hyperbolic  conservation  laws  is  the  inherent  presence 
of  discontinuities  (shocks)  in  the  solution  (Lebacque,  1996).  Moreover,  discontinuous  weak 
solutions  are  not  unique  for  hyperbolic  conservation  laws  and  entropy  conditions  must  be 
satisfied  to  obtain  physically  valid  solution  that  is  consistent  with  human  behavior  (such 
as  the  driver’s  ride  impulse)  (Ansorge,  1990;  Velan  and  Florian,  2002).  For  some  specific 
forms  of  the  equilibrium  flow-density  relationship,  the  LWR  model  can  be  solved  analytically 
(Haberman,  1977;  Whitham,  1974).  Recently,  the  analytical  solution  for  a  specific  class  of 
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LWR  model  was  derived,  which  assumed  that  the  flow-density  relationship  is  governed  by 
a  quadratic  function  throughout  the  density  regime  (Wong  and  Wong,  2002b).  However, 
from  many  observations,  this  assumption  is  too  strong  for  general  applications.  Therefore, 
it  is  natural  to  extend  the  method  to  the  case  with  a  more  general  flow-density  relationship. 
In  the  traffic  flow  literatnre,  it  is  not  nncommon  to  divide  the  range  of  the  density  into  a 
nnmber  of  regimes,  within  each  of  which  a  flow-density  cnrve  is  fitted  with  the  observed 
data  (Dick,  1966;  Edie,  1961;  May  and  Keller,  1968;  Underwood,  1961).  In  this  paper, 
we  assume  that  the  flow-density  relationship  is  concave  and  is  represented  by  a  piecewise 
qnadratic  fnnction,  with  any  two  adjacent  pieces  joining  continnously  at  a  critical  density  pc 
bnt  with  a  discontinuous  first  derivative  at  pc-  Thus  the  considered  flow-density  relationship 
is  only  Lipschitz  continnous  bnt  not  everywhere  differentiable.  For  such  type  of  flow-density 
relationships,  the  entropy  solution  still  exists  and  is  unique  for  general  bounded  variation 
initial  conditions  (Dafermos,  1972).  We  obtain  in  Section  2  explicit  formulas  for  the  entropy 
solntions  with  snch  flow-density  relationship  and  with  piecewise  linear  initial  condition  and 
piecewise  constant  bonndary  conditions.  In  Section  3  we  snmmarize  the  solution  procednre, 
concentrating  on  the  discussion  of  finding  the  earliest  time  when  the  waves  (characteristic 
lines  or  shocks)  from  the  previous  initial  condition  intersect  with  one  another  and  hence 
the  construction  of  the  entropy  solution  must  be  restarted  based  on  a  new  piecewise  linear 
initial  data.  In  Section  4  we  provide  nnmerical  examples  in  traffic  flows  to  demonstrate 
the  explicit  solntions  obtained  in  Section  2.  We  also  compare  these  explicit  solntions  with 
nnmerical  solntions  obtained  by  using  the  high  order  weighted  essentially  non-oscillatory 
(WENO)  schemes  (Jiang  and  Shu,  1996;  Zhang  et  ah,  2003).  Concluding  remarks  are  given 
in  Section  5. 
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2  Explicit  construction  of  the  entropy  solutions 


The  governing  equation  for  the  LWR  model  is  the  following  scalar  hyperbolic  conservation 
law 

Pt  +  q{p)x  =  0  (1) 

with  suitable  initial  and  boundary  conditions.  Here  p  G  (0,pmax)  is  the  density,  p^ax  is  the 
maximum  (jam)  density,  and  g(p)  is  the  traffic  ffow  on  a  homogeneous  highway,  which  is 
assumed  to  be  a  function  of  the  density  p  only  in  the  LWR  model.  More  specihcally,  the 
flow  g,  the  density  p  and  the  equilibrium  speed  u  are  related  by 

q{p)=u{p)p.  (2) 


In  this  paper,  the  flow  g(p)  is  considered  to  be  continuous,  piecewise  quadratic,  and  concave. 
Without  loss  of  generality,  we  will  concentrate  our  discussion  on  the  situation  where  the  flow 


q  is  dehned  by  two  different  quadratic  functions  in  different  regimes 


q{p)  = 


gi(p),  0<p<pc 
g2(p),  pc<p<pT, 


where 


Flux  I;  gi(p)  =  do  +  di  p  +  (i2P^;  Flux  II;  g2(p)  =  Cq  +  Ci  p  +  62  p^  (4) 


are  two  different  quadratic  functions,  which  are  continuous  at  the  junction  gi(pc)  =  g2(Pc), 
concave  in  each  piece  gi(p)  <  0  and  q'^^p)  <  0,  and  concave  also  at  the  junction  gj(pc)  > 
g2(Pc)-  A  typical  flow  in  this  setup  is  given  in  Figure  1.  The  general  situation  of  the  flow  q 
with  more  than  two  pieces  of  quadratic  functions  can  be  considered  with  the  same  recipe  to 
each  neighboring  pairs  of  quadratic  flow  functions. 

We  now  start  the  construction  of  explicit  solutions  to  the  conservation  law  (1)  with  such 
flows  g(p),  when  the  initial  condition  is  piecewise  linear.  We  will  hrst  ignore  the  boundary 
conditions,  and  will  leave  the  discussion  of  the  treatment  of  piecewise  constant  boundary 
conditions  to  Sections  2.3  and  2.4.  We  begin  with  the  generalized  Riemann  problem 


p(a;,0) 


ai  +  /5i  x,  a;  <  0 
0(2  +  P2X,  X  >  0 


(5) 
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Figure  1:  A  typical  flow  with  two  different  concave  quadratic  functions  joining  continuously 
at  a  critical  density  pc  with  discontinuous  and  decreasing  derivative  at  pc- 

We  also  assume  that,  for  the  x  range  we  are  considering,  the  initial  density  at  +  f3iX  is 
completely  contained  in  one  of  the  regimes  p  <  Pc  or  p  >  pc  for  i  =  1  and  2.  This  does 
not  lose  generality,  as  we  can  break  a  single  linear  function  into  two  pieces  as  in  (5)  when 
it  crosses  the  critical  density  pc-  We  also  remark  that  we  do  not  need  to  consider  the  case 
when  both  linear  functions  +  f3iX,  for  i  =  1,2,  are  contained  in  a  single  regime  p  <  Pc  or 
P  >  Pc  because  this  is  covered  by  the  results  in  (Wong  and  Wong,  2002b). 

Before  solving  the  generalized  Riemann  problem,  we  write  down  a  simple  but  important 
fact  that  we  will  heavily  use  in  the  sequel.  For  the  scalar  conservation  law  (1)  with  a  quadratic 
flux  q{p)  =  a  +  b p  +  c p^  and  a  linear  initial  condition  p{x,  0)  =  a  +  (3x,  the  solution  stays 
linear 

p{x,t)  =  a{t)  +  I3(t)  X  (6) 
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with 


a  -  bl3t  13 


l  +  2c/3t’  l  +  2cl3t' 

This  can  be  easily  obtained  by  the  method  of  characteristics  and  we  can  easily  check  that 
(6)  with  (7)  is  a  smooth  solution  to  the  conservation  law  (1)  with  the  initial  condition 
p(a;,  0)  =  a  +  {3  X,  until  it  becomes  singular  (i.e.  until  the  denominator  1  +  2c/5t  =  0).  This 
simple  fact  is  the  main  reason  that  enables  us  to  obtain  explicit  formulas  for  the  entropy 
solution.  The  solution  for  each  linear  piece  of  the  initial  condition  is  given  by  (6)-(7)  until 
neighboring  waves  interact  with  each  other. 

We  now  assume  that  the  x-axis  is  divided  into  a  number  of  intervals,  within  each  of 
them  the  initial  density  is  given  by  a  linear  function  p(a;,  0)  =  a -\-  j3x  which  is  completely 
contained  in  one  of  the  regimes  p  <  Pc  or  p  >  p^.  We  consider  the  solution  to  the  generalized 
Riemann  problem  (1)  with  the  initial  condition  (5)  for  each  of  the  inner  boundary  points 
separating  two  piecewise  linear  initial  conditions.  The  left  and  right  intervals  to  the  inner 
boundary  point  under  consideration  are  denoted  by 


e={xi,Xr),  e  =  {xi,Xr) 


respectively,  with  clearly  Xr  =  xi.  The  initial  condition  density  values  at  the  relevant  interval 
boundaries  are  denoted  by 

Pz  =  p(a;+,0),  pr=p{x;,Q)-  p,  =  p(x+,0),  p^  =  p(x7,0) 

see  Figure  2.  Notice  that  the  exact  solution  is  obtained  only  to  the  smallest  time  when  the 
waves  (characteristic  lines  or  shocks)  from  the  initial  condition  intersect  with  one  another.  At 
this  time  the  new  piecewise  linear  initial  condition  is  formed  and  the  procedure  is  repeated. 
We  consider  the  following  situations  separately. 

2.1  Propagation  of  a  shock 

If  Pr  Pc  Pi  and  the  intervals  e  and  e  belong  to  Flux  I  and  Flux  II  in  (4)  respectively, 
as  shown  in  Figure  2,  then  a  shock  satisfying  the  Lax  entropy  condition  (Lax,  1973)  is 
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Figure  2:  Propagation  of  a  right-moving  shock  from  a  nodal  point. 

generated  from  the  point  Xr  =  xi  and  will  move  to  the  right  or  left  along  a  curve  determined 
by  the  Rankine-Hugoniot  jump  condition.  If  the  shock  moves  to  the  right,  an  easy  way  to 
determine  the  location  of  the  shock  xi  -|-  Ax  after  time  At  is  through  conservation  in  the 
rectangular  region  with  (xi,  0)  and  (xi  +  Ax,  At)  as  the  end  points  of  a  diagonal  (see  Figure 
2).  Notice  that,  as  we  consider  only  the  time  At  smaller  than  the  smallest  time  when  the 
waves  (characteristic  lines  or  shocks)  from  the  initial  condition  intersect  with  one  another, 
we  may  safely  assume  that  the  left  and  top  boundaries  of  this  rectangle,  dQi  and  dQt,  belong 
to  Flux  I,  and  the  right  and  bottom  boundaries  of  this  rectangle,  dflr  and  dQb,  belong  to 
Flux  II. 

The  flux  at  the  left  boundary  dfli,  namely  the  number  of  vehicles  coming  from  the  left 
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boundary  into  the  region  hi  during  time  period  At  is 


pAt  j‘At 

fi=  qi\x=xrdt=  [do  +  di{ai{t)  +  Xr)  +  d2{ai{t)  +  Xrf]  dt  (8) 

Jo  Jo 

where  Q;i(t)  and  Pi(t)  are  given  by  (7)  with  a  and  f3  replaced  by  ai  and  (3i  in  the  initial 
condition  (5).  Likewise,  the  flux  at  the  right  boundary  dflr,  namely  the  number  of  vehicles 
leaving  the  right  boundary  from  the  region  hi  is 

pAt 

fr  /  Q2\x=xi+Ax  dt 


rAt 


eo  +  ei{a2{t)  + /32{t){xi  +  Ax))  +  e2{a2{t)  + /32{t){xi  +  Ax))"^]  dt.  (9) 


The  initial  number  of  vehicles  within  the  region  hi  at  time  t  =  0  is 

pxi-\-Ax 


h  = 


(0:2(0)  +  /92(0)a;)  dx 


fxi 


and  the  final  number  of  vehicles  within  the  region  at  time  t  =  At  is 

pXr-\-Ax 

ft=  (ai(At)  +  Pi{At)x)  dx. 

J  Xr 

From  the  flow  conservation  principle,  we  deduce  that 

fi  —  fr  +  fb  —  ft  =  0. 


(10) 


(11) 


(12) 


Using  the  explicit  formulas  (6)-(7),  we  obtain  from  (12)  the  explicit  equation  determining 
Ax  as 

Fi{At)  Ax‘^  +  F2{At)  Ax  +  Fo{At)  =  0  (13) 

where 

Fi(At)  =  {2At{d2  -  62)  (Pr-  -  Pl){A  -  Pi)  +  -  Xl){p^  -  Pi)  -  {Xr  -  Xi){pr  -  pi)}  / 

{2  [2d2At{pi  -  Pr)  FXi-  Xr]  [2e2At{'Pi  -  pj  FXi-  Xr]} 

F2{At)  =  {-{xi  -  Xr)[At{ei  +  2e2pr){pl  “  Pr)  -  {pi  -  Pr){xi  -  Xr)]  +  diAt{pi  -  Pr)[2e2At{pi 

+Xi  -  Xr]  -2d2At{pi  -  Pr)[Atei{pi  -  Pr)  -  pfxi  -  Xr)]}  / 

{[2d2At{pi  -  Pr)  +Xi-  Xr][2e2At{pi  -  Pr)  +Xi-  Xr]} 


FsiAt)  =  {At{{xi  -  Xr){At[el  +  4(do  -  eo)e2]ipi  -  Pr)  +  2(do  -  Cq  -  CiP;  -  e2p‘f){xi  -  Xr)  } 

+2d2{Af[el  +  4{do  -  eo)e2]{pl  -  pr){pl  -  Pr)  +  ‘^^t{e2[pipl{xi  -  Xr)  +  PrAi^r  -  Xl) 
-pKPI  -  Pr){xi  -  A)]  +  {do  -  eo  -  eipi){pi  -  Pr){xi  -  A)}  +  pl{xi  -  Xr){A  -  ^r)} 
-d\At{pi  -  pr)[2Ate2{A  -  Pr)  +  A-  Xr]  +2dipr{xi  -  Xr){2Ate2{pi  -  Pr)  +Xi-  Xr)}} 
/  {2[2d2At{pi  -  pr)  +Xi-  Xr][2e2At{pi  -  pj  +Xi-  Xr]} 


The  shock  trajectory  can  therefore  be  determined  by  solving  the  quadratic  equation  (13) 


Ax 


-F2{At)  +  VA 
2Fi(At) 


(14) 


where  A  =  F2{At)‘^  —  4Fi{At)F3{At).  If  the  root  Ax  determined  by  (14)  is  negative,  then 
the  assumption  of  a  right-moving  shock  is  incorrect.  In  this  case,  we  can  determine  the 
location  of  the  shock  xi  -|-  Ax  after  time  At  again  through  conservation  in  the  rectangular 
region  with  (Tz,0)  and  {xi  +  Ax,  At)  as  the  end  points  of  a  diagonal  (see  Figure  3).  We 
may  again  safely  assume  that  the  left  and  bottom  boundaries  of  this  rectangle,  dQi  and  dQb, 
belong  to  Flux  I,  and  the  right  and  top  boundaries  of  this  rectangle,  dQr  and  dQf,  belong 
to  Flux  II.  The  left  and  right  fluxes  fi  and  fr  are  changed  to 

rAt 

fl  /  dl\x=Xr  +  Ax  dt 

Jo 

At 

[do  +  di{ai{t)  +  I3i{t){xr  +  Ax))  +  d2{ai{t)  +  I3i{t){xr  +  Ax)Y]  dt  (15) 


and 


plAt  pAt 

fr=  q2\x=xidt=  [eo  +  ei{a2{t)  +  p2{t)A) +  e2{a2{t)  +  p2{t)A)^]  dt 

Jo  Jo 

respectively.  The  bottom  and  top  fluxes  are  changed  to 

PXr-\-Ax 


(16) 


A  = 


(q;i(0)  -|-  /3i(0)a;)  dx; 


(17) 


and 


pxi-\-Ax 

ft=  {a2{At)  +  (32{At)x)  dx. 

Jxi 


(18) 
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In  this  case,  the  flow  conservation  principle  (12)  yields  the  explicit  equation  determining  Ax 
as 

Fi(At)  +  F2{At)  Ax  +  Fs{At)  =  0  (19) 

where  we  can  easily  check  that  Fj(At)  =  —Fi{At),i  =  1,  2,  3.  Therefore  the  shock  trajectory 
can  be  determined  by  solving  the  quadratic  equation  (19) 


Ax  = 


-F^iAt)  -  y/A  -F^iAt)  +  ^/A 


(20) 


2Fi{At)  2Fi(At) 

where  A  =  F2{AtY  —  AFi{At)FY^t)  =  F2{AtY  —  4Fi(At)F3(At).  Notice  that  the  formulas 
(14)  and  (20)  are  identical  so  we  do  not  need  to  decide  a  priori  whether  the  shock  moves  to 
the  left  or  right. 


Figure  3:  Propagation  of  a  left-moving  shock  from  a  nodal  point. 
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2.2  Propagation  of  a  rarefaction  wave 

If  Pr  ^  Pc  ^  Pi  and  the  intervals  e  and  e  belong  to  Flux  II  and  Flux  I  in  (4)  respectively, 
then  a  rarefaction  wave  is  generated  and  three  new  intervals  ei,  62  and  es  are  created  at 
the  time  t  =  At,  as  shown  in  Figure  4.  We  again  consider  only  the  time  At  smaller  than 
the  smallest  time  when  the  waves  (characteristic  lines  or  shocks)  from  the  initial  condition 
intersect  with  one  another.  The  coordinates  of  the  four  nodes  serving  as  the  end  points  of 
the  three  intervals  ei,  62  and  63  at  time  At  can  be  determined  as 

Xi{At)  =  Xr  +  q2{Pr)At,  0:2  (At)  =  Xr  +  q2{Pc)At, 

X3{At)  =  Xr  +  q[{pc)At,  X4{At)  =  Xr  +  q[{pi)At. 

The  density  at  these  end  points  at  time  At  are  given  by 

p{xi{At),  At)  =  pr,  p{x2{At),  At)  =  pc,  p{xs{At),  At)  =  pc,  p{x4{At),  At)  =  pi 

and  the  density  is  linear  within  each  of  the  new  intervals  ej,  i  =  1,2,  3.  In  particular,  the 
density  within  e2  is  a  constant  p  =  pc- 

2.3  Boundary  conditions  from  the  highway  entrance 

In  this  subsection  we  discuss  the  boundary  condition  at  the  left  boundary  a;  =  0,  which  is 
the  highway  entrance. 

Looking  at  the  solution  (6)-(7),  we  can  see  that  the  general  solution  with  a  linear  initial 
condition  is  not  linear  in  t  for  hxed  x,  unless  /?  =  0,  in  which  case  the  solution  is  constant  in  t. 
Therefore,  within  our  piecewise  linear  (in  space)  framework,  we  can  only  consider  piecewise 
constant  boundary  conditions.  A  typical  boundary  condition  at  a;  =  0  is  shown  in  Figure  5. 

We  now  use  the  notation  in  the  previous  sections.  The  interface  is  at  Xr  =  xi  =  0,  and  the 
left  density  value  at  a;  =  0  is  given  by  the  boundary  condition.  Otherwise,  this  is  identical 
to  the  situation  studied  in  the  two  previous  subsections  for  an  internal  generalized  Riemann 
problem.  Again,  we  would  need  to  consider  the  situation  where  pr  and  the  linear  function  in 
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Figure  4:  Propagation  of  a  rarefaction  wave  from  a  nodal  point. 

the  first  interval  with  end  values  p;  and  belong  to  different  regimes  in  (4),  for  otherwise 
the  solution  is  the  one  obtained  in  (Wong  and  Wong,  2002b).  If  Pr  <  Pc  <  Pn  a  shock  is 
generated.  We  consider  only  the  situation  that  the  shock  moves  to  the  right.  Its  location 
Ax  after  time  At  is  given  by  (13),  except  that  the  left  flux  fi  given  by  (8)  is  simplihed  to 

fi  =  gi(pr)  At 

and  the  top  flux  ft  given  by  (11)  is  simplified  to 

ft  =  pr  Ax. 

Therefore,  the  coefficients  in  the  quadratic  equation  (13)  which  determines  the  shock  location 
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p 


Pi 


0 


P2 


t 


Figure  5:  A  typical  boundary  condition  at  the  entrance  a;  =  0. 

Ax  are  simplified  to 
Fi(At)  =  p^-pi 

F2{At)  =  2[ei{pi  -  +  2e2Pr{Pl  “  +  {pr  -  Pl){xi  -  Xr)] 

Fsi^At)  =  -At{[el+4:{do-eo)e2]{pi-Pr)At  +  2{do-eo-eiPi-e2p‘{)(xi-Xr) 

+  2diPr[2e2{Pi  -  Pr)^t  +  Xi-  Xr]  +  2d2pl[2e2{Pi  -  Pr)^t  +  Xi-Xr]}. 

When  Pr  A  Pc  A  Pi,  rarefaction  wave  is  formed.  The  formulas  are  the  same  as  those 
given  in  Section  2.2. 

2.4  Boundary  conditions  for  the  highway  exit 

In  this  subsection  we  discuss  the  boundary  condition  at  the  right  boundary  x  =  Xend,  which 
is  the  highway  exit. 
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We  consider  a  typical  exit  setup  with  a  traffic  signal,  which  alternates  between  green  and 
red  lights.  This  is  similar  to  the  piecewise  constant  boundary  condition  considered  for  the 
entrance  in  the  previous  subsection.  The  constant  values  of  the  density  p  for  the  green  and 
red  lights  are  P;  =  0  and  p;  =  Pmax,  respectively. 

When  p;  =  pmax,  corresponding  to  the  situation  of  a  red  light,  a  left-moving  shock  is 
formed.  We  again  only  consider  the  situation  Pr  <  Pc  <  P;,  corresponding  to  the  situation 
where  the  density  at  the  left  of  Xend  belongs  to  the  regime  of  Flux  I  in  (4).  The  value 
Ax  determining  the  shock  location  Xend  +  Ax  after  time  At  (recall  that  in  this  case  Ax  is 
negative)  is  given  by  (19),  except  that  the  right  flux  given  by  (16)  is  simplihed  to 

fr  =  0 

and  the  top  flux  ft  given  by  (18)  is  simplified  to 

ft  ~  Pmax  Ax. 

Therefore,  the  coefficients  in  the  quadratic  equation  (19)  which  determines  the  shock  location 
Xend  +  Ax  are  simplihed  to 

Fi{At)  =  pi-pr 

A(At)  =  -2[di{pi  -  Pr)At  +  2d2Pi{pi  -  Pr)At  +  (p;  -  Pr){xi  -  Xr)] 

Fz{At)  =  At{d\{pi  -  Pc) At  -h  2diPr{-xi  +  Xr)  -  2{-[eo  +Pz(ei  4-  e2Pi)]{xi  -  Xr) 

+do[2d2{pi  -  Pr)At  +  Xi  -  Xr]  +  ^2 [-2e2Pzpf  At  -4  2e2pfpr-At  -4  2eo(pr  -  Pi)At 
+  2e{pi{pr  -  pl)At  +  PrXl  -  plxr]}} 

When  p;  =  0,  corresponding  to  the  situation  of  a  green  light,  a  rarefaction  wave  is  formed. 
The  formulas  are  the  same  as  those  given  in  Section  2.2. 

3  Solution  procedure 

In  this  section  we  summarize  the  solution  procedure,  concentrating  on  the  discussion  of 
hnding  the  earliest  time  when  the  waves  (characteristic  lines  or  shocks)  from  the  previous 


14 


initial  condition  intersect  with  one  another  and  hence  the  construction  of  the  entropy  solution 
must  be  restarted  based  on  a  new  piecewise  linear  initial  data. 

Recall  our  assumption,  as  stated  in  Section  2,  that  the  x-axis  is  divided  into  a  number  of 
intervals,  within  each  of  them  the  initial  density  is  given  by  a  linear  function  p{x,  0)  =  a  +  j3  x 
which  is  completely  contained  in  one  of  the  regimes  p  <  Pc  or  p  >  pc-  We  consider  several 
cases  of  wave  interactions  separately  below,  and  then  take  the  smallest  time  from  these  cases, 
which  will  serve  as  the  time  to  restart  the  solution  procedure  with  a  new  piecewise  linear 
initial  data.  The  cases  being  considered  are:  natural  break  time  for  each  individual  linear 
piece,  and  the  interaction  between  waves  from  two  adjacent  nodes  (including  the  boundary 
nodes).  For  the  latter  case,  we  also  distinguish  between  the  situation  of  the  two  adjacent 
nodes  corresponding  to  a  shock  and  a  rarefaction  wave  respectively,  and  the  situation  of 
the  two  adjacent  nodes  both  corresponding  to  shocks.  Notice  that  the  situation  of  the  two 
adjacent  nodes  both  corresponding  to  rarefaction  waves  are  included  in  the  case  of  natural 
break  time  for  the  individual  linear  piece  between  the  two  nodes. 

Once  the  smallest  time  of  wave  interactions  is  found,  the  density  function  is  obtained 
at  this  time  using  the  formulas  in  Section  2  as  a  new  piecewise  linear  function,  and  the 
procedure  is  repeated,  until  the  hnal  desired  time  is  reached. 

3.1  Natural  break  time  for  each  individual  linear  piece 

For  an  element  e  in  which  the  initial  linear  density  prohle  is  an  increasing  function,  the 
natural  break  time,  see  Figure  6,  can  be  determined  by  the  following  formula  (Whitham, 
1974): 

Notice  that  under  our  assumption  (concave  quadratic  flux  q  and  increasing  linear  density 
p)  the  denomination  is  a  negative  constant,  hence  Te  is  a  positive  constant.  Indeed,  if  the 
interval  e  =  {xi,Xr)  and  the  initial  condition  density  values  at  the  interval  boundaries  are 
Pi  =  pixf  and  pr  =  p(a;7,0),  see  Figure  6,  and  assuming  that  the  flux  function  in  the 
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element  e  is: 


q{p)  =ao  +  aip  +  a2p^ 


then 


Te 


Xr  —  Xi  ^ 

2a2ipi-pr)  >  P^ 

oo  otherwise 


(21) 


Figure  6:  Natural  break  for  an  increasing  linear  prohle. 

3.2  A  shock  intersecting  with  its  adjacent  characteristic  line 

We  consider  the  situation  that  a  shock  from  the  element  boundary  Xr  =  xi  intersects  with 
the  characteristic  line  from  the  right  boundary  Xr  of  the  element  e,  see  Figure  7. 

Assuming  the  flux  function  in  the  element  e  is 

<?!  (p)  =  do  +  dl  p  +  (^2  P^ 
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Figure  7:  A  shock  intersecting  with  its  right  neighbor’s  characteristic  line 
and  the  flux  function  in  the  element  e  is 

q-iip)  =  eo  +  eip  +  62/, 

by  substituting  the  straight  line  function  of  the  characteristic  line  from  the  node  Xr  into 
the  shock  location  trajectory  (13)  for  the  shock  from  the  node  Xr  =  xi,  we  obtain  a  cubic 
equation  for  the  intersecting  time  t\ 

Tit^  +  +  Tsf  +  T4  =  0  (22) 

where,  with  (5i  =  Xr—xi,  S2  =  Xr  —  xi,  we  have 

Ti  =  — 2e2{pi  —  Pr)iPi  —  Pr)[^i  —  4:dod2  +  4(i2Co  +  Ci  +  4eie2p^  —  Ad2€,2P^  +  4e2P^  —  2di{ei  +  2e2p^)] 
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T2  =  {dl^  +  el^){pi  —  Pr)di  —  Aeie2{pi{pi  —  2p^)5i  +  Pr[~Pi{di  +  52)  +  +  <^2)]} 

2di{—ei{pi  —  Pr)5i  +  2e2{pi{pi  ~  ‘^Pr)^i  +  +  <^2)  +  Pr(2<^i  +  <^2)]}} 

+4{e2{  — ((io  —  Go){pi  —  Pr)52  +  ^2Pr{Pl{~‘^Pl  +  +  Pr[~Pr{‘55l  +  262)  +  Pr52] 

5-Pl[—Pr^2  +  2,pr{5i  +  52)]}}  +  d2{  —  {dQ  —  eQ){pi  —  Pr)5l  +  G2{Pl{.Pl  —  2^pI)5i  + 

Pr[—pf5i  —  PiPr52  +  Pr{2p^5i  +  Pr52)]}}} 

T3  =  2{(i2[~Pi<^l  (Pi  +  Pr)  +  Pr{Pl5i  +  Pr<^l  +  Pr<5l<52)]  ~  ~  60)^2  +  C?l[pZ<^l  ~  Pr(<^l  +  <^2)] 

ei[— pz^i  +  Pr{5i  +  52)]}  +  e2{pz5i(— Pz  +  3p^)  —  pf5i52  +  PiPr5i{5i  +  2^2) 

— p^[— 2p^5i(52  +  Pr(3(53  +  4(5i(52)]}} 

T4  =  [(p;  +  P^)(52  +  Pz^i  —  Pr.((5i  +  252)] 

We  can  then  determine  the  earliest  intersection  point  of  the  shock  and  the  characteristic  line 
by  finding  the  smallest  positive  root  of  eqnation  (22).  As  this  is  a  cnbic  eqnation  and  exact 
root  formulas  exist,  such  a  root  can  be  found  readily. 

A  symmetric  situation  is  when  a  shock  from  the  element  boundary  Xr  =  xi  intersects 
with  the  characteristic  line  from  the  left  boundary  xi  of  the  element  e,  see  Figure  8. 

The  procedure  to  obtain  the  intersecting  time  is  similar.  We  have  the  cubic  equation 

fit^  +  f2t^  +  Tsf  +  T4  =  0  (23) 

where,  with  =  Xr  —  xi,  62  =  Xr—xi,  we  have 

Ti  =  2d2[d‘l  —  2diei  +  +  4e2(do  ~  ^0)  +  4(i2Pz(di  —  ci)  +  4d2pf(d2  —  e2)](pz  —  Pr){pi  —  Pr) 

T2  =  -di(p;  -  pj(5i  -  [e\  +  4(cio  -  eo)e2](pz  -  pJ5i  -  4d2pi{-pf52  +  2p4p^5i  -  p;((5i  +  52)] 
+Pz[~3p^(5i  +  Pr52  +  P;(3(5i  +  252)]}  —  2(ii{— ei(p;  —  p,,)5i  +  2d2{Pr[Pr^l  ~  Pz('^l  +  *^2)] 
+Pz[~2p^(5i  +  Pi{25i  +  52)]}}  +  4d2{—{dQ  —  eo){pi  —  pr)52  +  e2[2p^(p;  —  p^)5i  +  piPi52 
+Pr(~PzPr<^l  +  PrPr^l  ~  Pl^2)]  +  ^l{Pr[Pr^l  ~  Pli^l  +  <^2)]  +  Pz[~2p^(5i  +  P;(2(5i  +  52)]}} 
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(Pi)  (Pr)(Pi)  (pj 


Figure  8:  A  shock  intersecting  with  its  left  neighbor’s  characteristic  line 

Tz  =  ‘^{^2{[—PlPr  +  PrPr  ^  +  Pr)]^l  ~  ^  ~  ^0)^2 

'^[Pr^l  ~  Pli^l  +  '^2)](ei  —  di)} 

+d2{—Pldi52  +  Pr[Prd‘i  +  Prdld2  —  P;(5i((5i  +  2^2)]  +  pl[—?>Pj.6\  +  Pi{^d\  +  45i52)]}} 
T4  =  5\[Pj.5i  +  {pi  +  Pr)d2  —  Pi{di  +  252)] 

3.3  Two  intersecting  adjacent  shocks 

We  now  consider  the  situation  that  two  adjacent  nodes  Xm  and  Xm+i  correspond  to  two 
shocks  which  will  intersect  at  time  At,  see  Figure  9.  This  is  the  most  difficult  case  since 
no  closed  form  formula  exists  for  the  intersecting  time  At.  We  will  therefore  resort  to  a 
nonlinear  equation  solver  such  as  the  Newton’s  method. 

We  assume  that  the  shock  location  of  the  shock  from  the  node  Xm  at  time  At  is  Xm  +  ^Xm 
and  denote  Gm{,^t)  =  Ax^.  Likewise,  the  shock  location  of  the  shock  from  the  node  Xm+i 
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Figure  9:  Two  intersecting  adjacent  shocks. 

at  time  At  is  Xm+i  +  Ax m+i  and  we  denote  Gm+i(At)  =  Ax^+i-  The  displacements  Gm{At) 
and  Gm+i{At)  are  determined  by  (13)  or  (19).  Therefore,  we  can  dehne  the  function 

S{At)  =  Xm+l  —  Xm  +  Gm+l{At)  —  Gm{At)  (24) 

which  measures  the  distance  between  the  two  shocks.  Clearly,  S'(O)  =  Xm+i  —  Xm  >  0,  and 
we  would  like  to  hnd  the  root  of  S{At),  which  corresponds  to  the  time  that  the  two  shocks 
intersect.  As  for  each  fixed  At,  S{At)  and  S' {At)  can  both  be  readily  computed,  we  can  easily 
set  up  a  Newton  iteration  to  solve  for  the  root  of  S{At)  with  At  =  0  as  the  initial  guess.  In 
the  Appendix,  we  will  show  the  uniqueness  of  the  intersecting  point  of  the  two  shock  curves 
before  the  natural  break  time  of  the  element  in  between,  if  they  do  intersect.  Therefore, 
S {At)  has  at  most  one  zero  before  the  natural  break  time  of  the  element  in  between,  which 
facilitates  the  convergence  of  the  Newton  iteration  procedure. 
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4  Numerical  examples 


In  this  section  we  provide  three  numerical  examples  to  illustrate  the  explicit  formulas  for  the 
entropy  solutions  obtained  in  the  previous  sections.  The  flow-density  relationship  is  given 
by  (3)-(4),  with 

(  -0.  V  +  lOOp,  0  <  p  <  50 

q{p)  =  <  -0.1p2  +  15p  +  3500,  50  <  p  <  100 

[  -0.024p2  -  5.2p  +  4760,  100  <  p  <  350 

see  Figure  10. 


Fundamental  diagram 


Figure  10;  Flow-density  relationship  used  in  the  numerical  examples. 


We  also  use  the  hfth  order  hnite  difference  WENO  scheme  (Jiang  and  Shu,  1996;  Zhang  et 
ah,  2003)  to  compute  the  solution  and  make  a  comparison.  The  purpose  of  this  comparison 
is  to  validate  the  computation  of  WENO  schemes  against  the  presumably  exact  entropy 
solutions  obtained  by  the  procedure  in  this  paper. 

Example  1.  Consider  a  homogeneous  highway  with  a  length  of  2  km.  Due  to  a  certain 
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incident  on  the  road,  an  initial  density  distribution  with  a  peak  of  150  veh/km  is  formed  as 
shown  in  Figure  11.  Assume  that  the  entrance  of  the  highway  is  blocked  and  no  traffic  is 
allowed  to  enter  from  the  upstream  end.  The  exact  solution  of  this  problem  can  be  worked 
out  using  the  procedure  in  this  paper,  shown  as  solid  lines  in  Figure  12,  in  comparison  with 
the  numerical  solution  obtained  by  the  WENO  scheme  using  N  =  200  uniform  grid  points, 
shown  as  circles.  We  can  see  that  the  two  results  agree  quite  well. 


Figure  11:  The  initial  density  prohle  for  Example  1. 

In  order  to  demonstrate  the  performance  of  the  solution  procedure  derived  in  this  paper, 
we  show  in  Figure  13  the  space-time  diagram  of  all  waves  (shocks  and  rarefaction  waves) 
involved  in  the  evolution  process  for  Example  1  until  t  =  3.0  min.  The  relevant  breaking 
points  and  density  values  at  nodes  are  also  shown  in  Table  1.  We  can  see  that  the  procedure 
works  well  for  this  example  with  relatively  few  restarts  to  reach  the  target  time  t  =  3.0  min. 
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Table  1:  Breaking  points  and  density  values  at  nodes  for  the  evolution  process  for  Example 
1  until  t  =  3.0  min. 


Time  (min) 

Element 

Xi  (km) 

Xr  (km) 

Pi  (veh/km) 

Pr  (veh/km) 

to  =  0 

1 

0.000 

0.167 

0.0 

50.0 

2 

0.167 

0.333 

50.0 

100.0 

3 

0.333 

0.500 

100.0 

150.0 

4 

0.500 

1.000 

150.0 

150.0 

5 

1.000 

1.167 

150.0 

100.0 

6 

1.167 

1.333 

100.0 

50.0 

7 

1.333 

1.500 

50.0 

0.0 

8 

1.500 

2.000 

0.0 

0.0 

ti  =  0.162 

1 

0.000 

0.271 

0.0 

0.0 

2 

0.271 

0.313 

82.4 

97.7 

3 

0.313 

0.466 

102.2 

150.0 

4 

0.466 

0.966 

150.0 

150.0 

5 

0.966 

1.140 

150.0 

100.0 

6 

1.140 

1.153 

100.0 

100.0 

7 

1.153 

1.347 

100.0 

50.0 

8 

1.347 

1.496 

50.0 

50.0 

9 

1.496 

1.771 

50.0 

0.0 

10 

1.771 

2.000 

0.0 

0.0 

t2  =  0.211 

1 

0.000 

0.307 

0.0 

0.0 

2 

0.307 

0.456 

102.9 

150.0 

3 

0.456 

0.956 

150.0 

150.0 

4 

0.956 

1.131 

150.0 

100.0 

5 

1.131 

1.149 

100.0 

100.0 

6 

1.149 

1.351 

100.0 

50.0 

7 

1.351 

1.545 

50.0 

50.0 

8 

1.545 

1.852 

50.0 

0.0 

9 

1.852 

2.000 

0.0 

0.0 

h  =  0.300 

1 

0.000 

0.358 

0.0 

0.0 

2 

0.358 

0.438 

124.0 

150.0 

3 

0.438 

0.938 

150.0 

150.0 

4 

0.938 

1.117 

150.0 

100.0 

5 

1.117 

1.142 

100.0 

100.0 

6 

1.142 

1.358 

100.0 

50.0 

7 

1.358 

1.633 

50.0 

50.0 

8 

1.633 

2.000 

50.0 

0.0 
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Table  1:  (continued) 


t4  = 

0.425 

1 

0.000 

0.412 

0.0 

0.0 

2 

0.412 

0.912 

150.0 

150.0 

3 

0.912 

1.096 

150.0 

100.0 

4 

1.096 

1.131 

100.0 

100.0 

5 

1.131 

1.369 

100.0 

50.0 

6 

1.369 

1.758 

50.0 

50.0 

7 

1.758 

2.000 

50.0 

23.1 

h  = 

0.667 

1 

0.000 

0.505 

0.0 

0.0 

2 

0.505 

0.862 

150.0 

150.0 

3 

0.862 

1.056 

150.0 

100.0 

4 

1.056 

1.111 

100.0 

100.0 

5 

1.111 

1.389 

100.0 

50.0 

6 

1.389 

2.000 

50.0 

50.0 

te  = 

1.274 

1 

0.000 

0.737 

0.0 

0.0 

2 

0.737 

0.954 

150.0 

100.0 

3 

0.954 

1.061 

100.0 

100.0 

4 

1.061 

1.439 

100.0 

50.0 

5 

1.439 

2.000 

50.0 

50.0 

h  = 

1.600 

1 

0.000 

0.900 

0.0 

0.0 

2 

0.900 

1.033 

100.0 

100.0 

3 

1.033 

1.467 

100.0 

50.0 

4 

1.467 

2.000 

50.0 

50.0 

ts  = 

1.778 

1 

0.000 

1.019 

0.0 

0.0 

2 

1.019 

1.481 

100.0 

50.0 

3 

1.481 

2.000 

50.0 

50.0 

tg  = 

2.333 

1 

0.000 

1.528 

0.0 

0.0 

2 

1.528 

2.000 

50.0 

50.0 

tio  = 

2.711 

1 

0.000 

2.000 

0.0 

0.0 

til  = 

3.000 

1 

0.000 

2.000 

0.0 

0.0 
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Example  2.  We  consider  a  long  homogeneous  freeway  of  length  20  km.  The  entrance 
density  is  50  veh/km.  Due  to  an  incident  near  the  downstream  end  of  the  freeway,  the  traffic 
density  profile  shown  in  Figure  14  is  formed  in  which  a  jam-packed  condition  of  5  km  long 
occurs  from  10  to  15  km  measured  from  the  upstream  entrance  of  the  freeway.  In  order  to 
release  the  traffic  jam  condition  downstream,  the  authority  blocks  the  freeway  entrance  for 
10  min,  after  which  traffic  is  released  again  from  the  freeway  entrance  at  the  capacity  density 
of  75  veh/km.  After  20  min,  the  entrance  flow  returns  back  to  normal  with  a  density  of  50 
veh/h.  The  variation  of  traffic  density  at  the  upstream  entrance  of  the  freeway  is  illustrated 
in  Figure  15.  The  period  of  analysis  is  2  hours. 

The  exact  solution  of  this  problem  is  worked  out  using  the  procedure  in  this  paper,  shown 
as  solid  lines  in  Figure  16,  in  comparison  with  the  numerical  solution  obtained  by  the  WENO 
scheme  using  N  =  200  uniform  grid  points,  shown  as  circles.  We  can  again  see  that  the  two 
results  agree  quite  well. 

Again,  in  order  to  demonstrate  the  performance  of  the  solution  procedure  derived  in  this 
paper,  we  show  in  Figure  17  the  space-time  diagram  of  all  waves  (shocks  and  rarefaction 
waves)  involved  in  the  evolution  process  for  Example  2  until  t  =  120  min.  The  relevant 
breaking  points  and  density  values  at  nodes  are  also  shown  in  Table  2.  We  can  see  that  the 
procedure  works  well  for  this  example  with  relatively  few  restarts  to  reach  the  target  time 
t  =  120  min. 

Example  3.  This  example  has  the  same  initial  condition  and  entrance  boundary  condition 
as  those  in  Example  2.  However,  at  the  exit  boundary,  a  traffic  signal  is  installed,  with  a 
repeated  pattern  of  2  minutes  green  light  followed  by  1  minute  red  light. 

As  before,  the  exact  solution  of  this  problem  is  worked  out  using  the  procedure  in  this 
paper,  shown  as  solid  lines  in  Figure  18,  in  comparison  with  the  numerical  solution  obtained 
by  the  WENO  scheme  using  N  =  400  uniform  grid  points,  shown  as  circles.  We  see  once 
again  that  the  two  results  agree  quite  well. 
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Table  2:  Breaking  points  and  density  values  at  nodes  for  the  evolution  process  for  Example 
2  until  t  =  120  min. 


Time  (min) 

Element 

Xi  (km) 

Xr  (km) 

Pi  (veh/km) 

Pr  (veh/km) 

to  =  0 

1 

0.000 

10.000 

50.0 

50.0 

2 

10.000 

15.000 

350.0 

350.0 

3 

15.000 

18.571 

350.0 

100.0 

4 

18.571 

19.286 

100.0 

50.0 

5 

19.286 

20.000 

50.0 

0.0 

ti  =  0.714 

1 

0.000 

0.952 

0.0 

0.0 

2 

0.952 

9.841 

50.0 

50.0 

3 

9.841 

14.738 

350.0 

350.0 

4 

14.738 

18.452 

350.0 

100.0 

5 

18.452 

18.512 

100.0 

100.0 

6 

18.512 

19.345 

100.0 

50.0 

7 

19.345 

20.000 

50.0 

50.0 

t2  =  6.429 

1 

0.000 

8.571 

0.0 

0.0 

2 

8.571 

12.643 

350.0 

350.0 

3 

12.643 

17.500 

350.0 

100.0 

4 

17.500 

18.036 

100.0 

100.0 

5 

18.036 

19.821 

100.0 

50.0 

6 

19.821 

20.000 

50.0 

50.0 

ts  =  8.571 

1 

0.000 

8.571 

0.0 

0.0 

2 

8.571 

11.857 

350.0 

350.0 

3 

11.857 

17.143 

350.0 

100.0 

4 

17.143 

17.857 

100.0 

100.0 

5 

17.857 

20.000 

100.0 

50.0 

o 

o 

o 

o 

1 

0.000 

8.571 

0.0 

0.0 

2 

8.571 

11.333 

350.0 

350.0 

3 

11.333 

16.905 

350.0 

100.0 

4 

16.905 

17.738 

100.0 

100.0 

5 

17.738 

20.000 

100.0 

52.5 

t5  =  15.143 

1 

0.00 

0.429 

75.0 

50.0 

2 

0.429 

5.143 

50.0 

50.0 

3 

5.143 

8.571 

50.0 

0.0 

4 

8.571 

9.448 

350.0 

350.0 

5 

9.448 

16.048 

350.0 

100.0 

6 

16.048 

17.310 

100.0 

100.0 

7 

17.310 

20.000 

100.0 

58.5 
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Table  2:  (continued) 


^6  — 

18.182 

1 

0.000 

0.682 

75.0 

50.0 

2 

0.682 

8.182 

50.0 

50.0 

3 

8.182 

8.333 

350.0 

350.0 

4 

8.333 

15.541 

350.0 

100.0 

5 

15.541 

17.056 

100.0 

100.0 

6 

17.056 

20.000 

100.0 

60.7 

h  = 

19.231 

1 

0.000 

0.769 

75.0 

50.0 

2 

0.769 

7.949 

50.0 

50.0 

3 

7.949 

15.366 

350.0 

100.0 

4 

15.366 

16.969 

100.0 

100.0 

5 

16.969 

20.000 

100.0 

61.3 

ts  = 

30.000 

1 

0.000 

1.667 

75.0 

50.0 

2 

1.667 

5.678 

50.0 

50.0 

3 

5.678 

13.571 

306.2 

100.0 

4 

13.571 

16.071 

100.0 

100.0 

5 

16.071 

20.000 

100.0 

65.6 

tg  = 

44.742 

1 

0.000 

0.698 

50.0 

50.0 

2 

0.698 

2.895 

69.0 

50.0 

3 

2.895 

11.115 

264.1 

100.0 

4 

11.115 

14.843 

100.0 

100.0 

5 

14.843 

20.000 

100.0 

68.4 

tio  = 

:  54.000 

1 

0.000 

1.195 

50.0 

50.0 

2 

1.195 

9.571 

245.7 

100.0 

3 

9.571 

14.071 

100.0 

100.0 

4 

14.071 

20.000 

100.0 

69.5 

til  = 

:  61.319 

1 

0.000 

8.352 

231.9 

100.0 

2 

8.352 

13.462 

100.0 

100.0 

3 

13.462 

20.000 

100.0 

70.1 

tl2  = 

111.429 

1 

0.000 

9.286 

100.0 

100.0 

2 

9.286 

20.000 

100.0 

72.2 

^13  = 

120.000 

1 

0.000 

8.571 

100.0 

100.0 

2 

8.571 

20.000 

100.0 

72.4 
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Table  3:  CPU  time  comparison  between  the  solution  procedure  in  this  paper  and  the  hfth 
order  WENO  scheme. 


Example 

Analytical  (sec) 

WENO  (sec) 

N 

ending  time  t  (min) 

Example  1 

1.669x10“^ 

0.3484 

200 

3 

Example  2 

2.675x10“^ 

0.4609 

200 

120 

Example  3 

7.933x10-^ 

2.045 

400 

120 

Finally,  in  Table  3  we  compare  the  CPU  time  used  by  the  procedure  in  this  paper  and  by 
the  hfth  order  WENO  scheme  with  an  appropriate  number  of  uniform  grid  points  to  achieve 
a  comparable  accuracy,  for  the  three  examples  shown  above.  We  can  see  clearly  that  the 
exact  solution  procedure  developed  in  this  paper  is  much  less  CPU  time  intensive  than  the 
hnite  difference  WENO  method. 

5  Conclusions 

We  have  developed  in  this  paper  a  procedure  to  explicitly  compute  the  entropy  solution  for  a 
Lighthill-Whitham-Richards  traffic  how  model  with  with  a  how-density  relationship  which  is 
piecewise  quadratic,  continuous,  concave,  but  not  diherentiable  at  the  junction  points,  and 
with  piecewise  linear  initial  condition  and  piecewise  constant  boundary  conditions.  The  pro¬ 
cedure  usually  involves  very  few  restarts  and  hence  rapid  computation  for  realistic  traffic  how 
situations,  using  very  small  CPU  time  comparing  with  a  high  order  hnite  diherence  WENO 
scheme.  Three  representative  numerical  examples,  involving  various  initial  and  boundary 
conditions  including  that  for  traffic  control  signals,  are  used  to  demonstrate  the  efficiency 
and  ehectiveness  of  this  solution  procedure.  The  solutions  obtained  agree  very  well  with 
those  obtained  by  the  hnite  diherence  WENO  scheme  which  uses  much  more  CPU  time. 
In  future  work  we  will  attempt  to  remove  the  condition  of  continuity  for  the  how-density 
relationship  in  order  to  better  model  certain  realistic  traffic  hows  such  as  the  famous  A-waves. 
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A  Appendix 

In  this  Appendix  we  study  in  more  detail  the  situation  discussed  in  Section  3.3,  that  two 
adjacent  nodes  correspond  to  two  shocks  which  will  intersect.  Referring  to  Figure  19,  we 
assume  that  two  shocks  Si(t)  and  5*2 (t),  starting  from  the  two  adjacent  nodes  xi  and  Xr, 
intersect  hrst  at  time  t*.  The  inner  element  e  has  two  edge  densities  p;  and  at  the  foot  of  the 
shocks  5'i(0)  and  5*2  (0).  We  further  assume  that  the  intersecting  time  t*  is  before  the  natural 
break  time  of  the  element  e  given  by  (21);  otherwise,  we  do  not  need  to  be  concerned  with 
the  intersection  time  of  the  two  shocks  anyway  because  the  solution  procedure  has  already 
been  restarted  before  the  two  shocks  have  the  chance  to  intersect. 

The  shock  curves  Si(t)  and  5*2 (^),  considered  individually,  may  certainly  be  dehned  be¬ 
yond  t  =  t*.  As  our  procedure  to  determine  the  intersection  time  t  =  t*  in  Section  3.3  is 
algebraic,  using  a  Newton  iterative  solver  to  hnd  the  root  of  the  nonlinear  equation  (24),  we 
would  like  to  prove  that  this  nonlinear  equation  has  at  most  one  root  before  the  natural  break 
time  of  the  element  e,  in  order  to  facilitate  the  convergence  of  the  Newton  iteration  proce¬ 
dure.  The  uniqueness  of  the  root  of  (24)  is  equivalent  to  the  uniqueness  of  the  intersecting 
point  t  =  t*  oi  Si{t)  and  5*2 (^)- 

We  prove  this  uniqueness  by  showing  the  monotonicity  of  the  density  values  pi(t)  (the 
right  value  of  the  left  shock  Si(t))  and  P2(t)  (the  left  value  of  the  right  shock  S2(t)),  see 
again  Figure  19.  By  continuity  of  the  density  between  the  two  adjacent  shocks,  we  clearly 
have  pi{t*)  =  P2{t*)-  Furthermore,  as  the  density  prohle  within  the  element  e  is  linear,  and 
the  characteristic  lines  between  the  two  adjacent  shocks  do  not  intersect,  it  is  easy  to  see 
geometrically  that  pi(t)  and  P2(t)  must  be  monotone  with  opposite  trends  (that  is,  either 
Pi{t)  is  monotonically  increasing  and  P2{t)  is  monotonically  decreasing,  or  vice  versa).  We 
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will  also  prove  this  fact  algebraically  below.  With  this  fact,  we  have  clearly  established  the 
uniqueness  of  the  intersecting  point  t  =  t*  of  the  two  shock  curves  Si(t)  and  5*2 (t). 

Assume  that  the  governing  flux  of  the  element  e  is 

q{p)  =  e2p^  +  eip  +  cq 


we  have,  by  (6)-(7),  the  formulas  for  pi(t)  and  P2(t): 

PlXr  -  prXl  +  ei{pi  -Pr)t-  {pi 


plit)  = 

P2it)  = 


Xr  -Xl  -  262  (p;  “ 
PlXr  -  prXl  +  ei{pi  -  p^)t  -  {pi 


Xr-Xi-  262  (Pi 


-Pr)iSl(t)  +Xi) 
Pr)t 

-pj(5'2(t)  +Xr) 
Pr)t 


Therefore  we  have 


,  ^  {Pl  -  Pr){(ei  +  2e2P;)(Tf  -  Xr)  +  262 (P;  "  Pr)Sl{t)  +  [262 (p^  "  Pl)t  +  W  -  Xi]S[{t)} 

[262 (Pz  -pr)t  +  Xi-  Xr]‘^ 

^  {Pl  -  pr){{ei  +  2e2'pr){Xi  -  Xr)  +  262 (p;  "  Pr)S2{t)  +  [262 (p^  "  Pl)t  +  W  - 

[262 (Pz  -'pr)t  +  Xi-  Xrf 

On  the  one  hand  S[(t)  >  q'{pi(t))  and  5*2 (f)  <  g'(p2(t)),  by  the  entropy  condition.  That 
is  to  say,  S[(t)  >  262Pi(t)  +  6i  and  5'2(t)  <  262P2(t)  +  Ci. 

On  the  other  hand  [262(Pr  ~  Pi)'^  +  Xr  —  Xi]  >0  since  we  assume  that  time  t  is  less  than 
the  natural  break  time  of  the  element  given  by  (21). 

We  now  have 


(6i  +  2e2Pi){xi  -  Xr)  +  262(Pz  “  Pr)Sl{t)  +  [262(p^  -  Pi)t  +  Xr  -  Xl]S[{t) 

>  (6i  +2e2Pi){xi  -  Xr)  +  262(Pz  +  [262 (p^  -pi)t  +  Xr  -  T;]  (262Pl  (t)  +  6i) 


=  (6i  +  2e2Pi){xi  -  Xr)  +  262 (Pz  “  Pr)Sl{t) 


+  [262(p^  -  Pi)t  +  Xr-  Xl]  (  262 


PlXr  -  PrXl  +  6i(pz  -  p^)t  -  {pi  -  pr){Si{t)  +  Xj) 
Xr  -  Xl  -  262(Pz  -Pr)^ 


+  ei 


(6i  +  262Pz)(Tz  -Xr)+  262 (p;  -  pJS'i(t) 

+262  [PzW  -  Pr+  +  ei(pz  -  pjt  -  (Pz  -  p^)(S'i(f)  +Xi)]  +  [262  (p^  -  Pi)t  +  Xr  -Xi\ei 
0 
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and 


(ei  +  2e2Pr){xi  -  Xr)  +  2e2(Pi  -  'Pr)S2{t)  +  [2e2(Pr  -  'PlY  +  Xr-  Xl]S'2{t) 

<  (ei  +  2e2Pr){xi  -  Xr)  +  2e2{pi  -  Pr)S2(t)  +  [2e2{Pr  -  Pl)t  +  Xr-  Xl]{2e2P2(t)  +  ei) 

=  (ei  +  2e2Pr){xi  -  Xr)  +  2e2{Pi  “  A)S2(t) 

+  I2e,(p,-  -  +  M 


V 


Xr  -Xi  -2e2{pi  -  Pr)t 


+  ei 


=  (ei  +  2e2Pr){xi  -  Xr)  +  2e2{Pi  “  Pr)S2(t) 

+  [2e2{Pr  -  Pl)t  +  Xr-  Xl]ei  +  2e2[PiXr  -  PrXi  +  ei(p;  -  p^)t  -  {Pi  -  p^){S2(t)  +  Xr)] 


=  0 


Therefore  p'i{t)  <  0,  P2{t)  >  0  if  pi  >  pp,  and  p[{t)  >  0,  P2{t)  <  0  if  pi  <  Pr-  The  desired 
monotonicity  is  proven. 
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Figure  12:  The  exact  entropy  solution  obtained  by  the  procedure  in  this  paper  (solid  line) 
and  the  numerical  solution  obtained  by  using  a  hfth  order  WENO  finite  difference  scheme 
with  N  =  200  uniform  grid  points  (circles).  Example  1.  Top  left:  t  =  0.2112  min;  top  right: 
t  =  0.4245  min;  bottom  left:  t  =  1.060  min;  bottom  right:  t  =  1.600  min. 
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Figure  13:  The  time-space  diagram  for  all  waves  (shocks  and  rarefaction  waves)  involved  in 
the  evolution  process  for  Example  1  until  t  =  3.0  min. 
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Figure  14:  The  initial  density  profile  for  Example  2. 
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Figure  15:  The  density  variation  at  the  entrance  in  time  for  Example  2. 
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Figure  16:  The  exact  entropy  solution  obtained  by  the  procedure  in  this  paper  (solid  line) 
and  the  numerical  solution  obtained  by  using  a  fifth  order  WENO  finite  difference  scheme 
with  N  =  200  uniform  grid  points  (circles).  Example  2.  Top  left:  t  =  10  min;  top  right: 
t  =  30  min;  bottom  left:  f  =  60  min;  bottom  right:  f  =  90  min. 
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Figure  17:  The  time-space  diagram  for  all  waves  (shocks  and  rarefaction  waves)  involved  in 
the  evolution  process  for  Example  2  until  t  =  60  min  (top)  and  until  t  =  120  min  (bottom). 
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Figure  18:  The  exact  entropy  solution  obtained  by  the  procedure  in  this  paper  (solid  line) 
and  the  numerical  solution  obtained  by  using  a  hfth  order  WENO  hnite  difference  scheme 
with  N  =  400  uniform  grid  points  (circles).  Example  3.  Top  left:  t  =  3  min;  top  right: 
t  =  30  min;  bottom  left:  t  =  60  min;  bottom  right:  t  =  90  min. 


42 


43 


